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Abstract. An efficient method for calculating inclusive conventional and prompt atmospheric leptons fluxes 
is presented. The coupled cascade equations are solved numerically by formulating them as matrix equation. 
The presented approach is very flexible and allows the use of different hadronic interaction models, realistic 
parametrizations of the primary cosmic-ray flux and the Earth’s atmosphere, and a detailed treatment of particle 
interactions and decays. The power of the developed method is illustrated by calculating lepton flux predictions 
for a number of different scenarios. 


1 Introduction 

Cosmic rays entering the Earth’s atmosphere produce a 
multitude of secondary particles in interactions with air 
nuclei. Some of the secondary particles decay into muons 
and neutrinos, which are not absorbed in the atmosphere 
and can reach particle detectors at ground level. The spec¬ 
tra of these leptons contains not only information about the 
primary cosmic rays, but also about the particle physics of 
their production and the properties of the traversed atmo¬ 
sphere. Furthermore, searches for high-energy neutrinos 
from astrophysical sources have to cope with a large flux 
of atmospheric leptons as background. A better under¬ 
standing of this flux, in particular its dependence on zenith 
angle and the changing properties of the atmosphere will 
help to develop improved methods to identify astrophys¬ 
ical neutrino fluxes and also contribute to a better under¬ 
standing of hadronic interactions at high energy. 

Many calculations of atmospheric lepton fluxes have 
been carried out since the early 1960’s (e.g. H|, see 12 for 
a review). However, most of them incorporated approx¬ 
imations in solving the cascade equations that lead to in¬ 
creased uncertainties on the relation between the predicted 
fluxes and the physical input parameters, or the calcula¬ 
tions could only be carried out in detail for just one or a 
few parameter/model combinations due to the large CPU 
time requirements. In this work we reduce the uncertain¬ 
ties related to the calculation method to a minimum by 
developing a numerical method with a level of detail com¬ 
parable with Monte Carlo calculations 00. In addition, 
the contributions of heavy flavor mesons and resonances to 
the flux of atmospheric leptons is accounted for in detail. 
The code is made publicly availably] 


a e-mail: anatoli.fedynitch@cem.ch 

* https ://github. com/afedy nitch/MCEq 


Our approach is based on the numerical solution of the 
coupled cascade equations, which have been rewritten into 
a matrix form to make use of modern implementations of 
linear algebra algorithms. While providing superior preci¬ 
sion at very high energies where Monte Carlo methods are 
often statistically inefficient, the high performance of the 
algorithm allows us to perform calculations for many in¬ 
put parameter and model assumptions in a very short time. 
On an average portable computer it takes a few seconds to 
calculate lepton fluxes, while keeping most of relevant pa¬ 
rameters accessible for users and easy to modify according 
to the current application. 

2 Coupled cascade equation 

The cascade equations for particle h can be written for one 
discrete energy bin E, 
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It is part of a system of coupled ordinary differential equa¬ 
tions, describing the evolution of the flux <f> of particles as 
a function of the atmospheric slant depth 
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For an observation height ho, X(ho) is computed along the 
trajectory / of the cascade core through the atmosphere us- 
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ing the mass density p, which is typically a function of 
the atmospheric height. The behavior of the particle cas¬ 
cade is driven by the competition of two source terms ©’ 
( |Td| and two sink terms ( fla| , ( fib} . The interaction length 
= »i air /o-™ l air (Ei ) in units g/cm 2 @ is independent 
of the slant depth and only varies slowly with energy due 
to the inelastic particle-air cross-section. The decay length 
A h , E (X) = cThEjp a j r (X)/mi, is proportional to the life¬ 
time Th of particle h and can vary by orders of magnitude 
due to the relativistic time dilation. In our approximation 
new particles are created along the shower trajectory in 
hadronic interactions in i fi~cj ) or decays in ( fld| ). with en¬ 
ergy conservation restricting the range of possible source 
particles to energies £/. equal or greater than £,. 

The interactions coefficients of particle / producing 
particle h, ci(E k )-+h(E t ), are obtained from hadronic interac¬ 
tions models by histogramming the particle yield as func¬ 
tion of x^ h = Ej/Ek for l -air collisions. Suitable interac¬ 
tion models include, for example, QGSJET-II-04 (7) and 
EPOS LHC {8j, and the upcoming versions of SIBYLL 
|9l and DPMJET-III itTOl . Alternatively one can directly 
extrapolate results of fixed-target experiments on light nu¬ 
clear targets using scaling arguments. 

While in some cases it is possible to find analytical 
expressions for the decay coefficients di(E k )-*h(Ei), see 0, 
numerical simulations have to be used for describing com¬ 
plex decays accurately. We have tabulated the decay coef¬ 
ficients as function of E /. based on simulations with the 
Monte Carlo event generator Pythia 8 HUH3 that in¬ 
cludes also rare decay channels and accounts for the effect 
of electroweak matrix elements where applicable. 

The full system of equations of the hadronic cascade 
can be obtained by writing Eq. (1) for all possible types of 
hadrons and leptons. In the following we will concentrate 
on high lepton energies and neglect the interaction and/or 
decay terms of neutrinos and muons. To obtain, for exam¬ 
ple, the flux of atmospheric neutrinos at the surface, one 
needs to solve the full system taking into account the non¬ 
linear X dependence of A dec and the non-analytic forms of 
the particle spectra serving as input for the interaction co¬ 
efficients 


10 1() GeV. We explicitly include more than 50 types of 
mesons, baryons and leptons. Together with some addi¬ 
tional technical groups the dimension of <l> is then ~ 6000. 

The reciprocal coefficients 1 /A of interaction and de¬ 
cay lengths are arranged in diagonal matrices 
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The decay length matrix A* c is constructed analogously 
usin g * h deC ' El = 4 fCrE .(Y)/p air (Y) to factorize out the de¬ 
pendence on the air density. Sub-matrices containing the 
interaction coefficients are defined as 
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and sub-matrices for the decay D^i, are constmcted in a 
similar way. The full interaction and decay matrices C and 
D are built from these sub-matrices according to the order 
of particle types in <I> 
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Using the definitions above, the matrix form of the coupled 
cascade equations can be written as 



(-1 + C)K m , + —1 + D)A dl 

p(X) 
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2.2 Short-lived particles 


2.1 Matrix form 


An efficient numerical computing scheme can be found 
by rewriting the cascade equations into matrix form. We 
group the different 0*(£,) into a column vector <J) by writ¬ 
ing blocks for each particle type for the discrete energy 
spectra 

®p(E i) 


0> = 


®p(En) 

<K(E 0 ) 


(3) 


The energy grid E, = 50GeV • 10 ,/iv is logarithmically 
spaced with roughly 8 bins per decade of energy across 
the energy range of the calculation between 50GeV and 


One goal of this work is to accurately take into account 
contributions of heavy flavor mesons and resonances to 
the flux of atmospheric leptons. Their short decay lengths 
introduce quickly decaying modes in Eqs. ([lb| and ( |Tdj >, 
which appear as eigenvalues of the decay matrix D with a 
large modulus of the negative real part. Fig. [T] illustrates 
the large difference between the decay lengths of conven¬ 
tional and charmed mesons. The energy dependence origi¬ 
nates from time dilation. If only pions and kaons would be 
considered as intermediate mesons, the cascade equations 
would become interaction dominated above 1 TeV. The de¬ 
cay would be a slow process and the equations could be 
easily integrated using a moderate step size of <9(1 g/cm 2 ). 
If short-lived particles are included, the choice of the ap¬ 
propriate step size is driven by the smallest eigenvalue, 
avoiding oscillations of fluxes around zero that are un¬ 
physical (<T > 0). For this reasons the equation system 
becomes a stiff numerical problem involving step sizes of 
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Figure 1. Decay lengths A dec for a subset of hadrons, evaluated 
at ham = 8 km. Superimposed is the interaction length A int of K ± . 


<9(1(T 5 g/cm 2 ). From a performance point of view this ap¬ 
proach is unreasonable since the integration can run up to 
values of X ~ <9(10 4 ) g/cm 2 . 

To reduce the stiffness of the equation system, we in¬ 
troduce the resonance approximation. In the resonance ap¬ 
proximation very short-lived particles (resonances), e.g. // 
or p mesons, decay immediately after their creation at the 
vertex. For each particle h this approximation is valid in a 
regime, where 


A 


h 

dec 


« A? 


(9) 


The coefficients for the chained production of a long lived 
secondary particle l via production and immediate decay 
of the resonance //, neglecting its interactions, can be writ¬ 
ten as 

= -O17—»/7 ' C/2—>77- (10) 

To understand this result, let i]"" be the vector containing 
the fluxes of all resonances, which are created during the 
integration step n. Using the matrix notation and forward 
Euler integration we can write 


rC = C^int <b • AX„. (11) 

C™ denotes a (k x c/®) matrix similar to C defined in 
Eq. (|7j, which contains production coefficients for reso¬ 
nances in interactions of hadrons. According to the ap¬ 
proximation, all created resonances have to decay into or¬ 
dinary particles within the same integration step. By writ¬ 
ing this condition for a single resonance type k 

Uk,n+ 1 = 0 — Tjf, n ~dlk dk.n ' AX n , 

A dec,eff 

we obtain the effective decay length 

A lc,eff = AX n=p(X)Al ceff . ( 12 ) 


Using APj[ ece ff instead of the true decay length for short¬ 
lived resonances we make sure that all particles decay af¬ 
ter one integration step in X without having to change the 
numerical treatment of the cascade equations. It should be 
noted that A^ ece ff does not depend on the properties of the 


resonance. The contribution to the ordinary particle flux 
due to resonance decay is then 


A®^=^A Z, eff ri- 

= D™ h r,. 


AX n 
PQ 0 


(13) 


where D r *l >h is a (of® x k) matrix, containing decay coeffi¬ 
cients of resonances into hadrons and leptons. By inserting 
Eq.© in m we replicate the expression from Eq. ( fT0| > 
for the production of particles via intermediate resonances 


A®-"r = (D res , 

n +1 v 77 —>h 

— R Kin, 


■ CZ n )\n, ■ AX„ 
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and define the square (d,,, X d,,,) resonance matrix R. 
Chained decays proceeding through two or more reso¬ 
nances are governed by additional left multiplications of 
decay matrices. Extending the matrix form of the cascade 
equations ([8]i with intermediate resonance production re¬ 
sults in 


=(-l + C + R)A inl <t> 
dX 

+ -^r(-l + D)A dec <h. (15) 

P(X) 

At high energies, where 

Adec ~ Aj n t, ( 16 ) 

the interaction of resonances becomes important. We in¬ 
troduce the parameter t m i X = Adec(E)/Aj nt {E), which is a 
threshold value, separating the energy regime where the 
particle can be treated as resonance from a regime where 
it has to be a full member of the cascade and listed in <D. 
A reasonable value is t mix = 0.05. This complicates some¬ 
what the procedure to fill the C, D and R matrices, where 
for each particle the individual threshold has to be taken 
into account as cut in row and/or column. 


3 Calculation input 

3.1 Initial state 

The initial state of the cascade equation is the flux of cos¬ 
mic rays at the top of the atmosphere. Since the calcula¬ 
tion relies on the properties of the average air shower, the 
superposition theorem is sufficient to model the flux and 
composition. A cosmic ray nucleus is modeled as Z pro¬ 
tons and A - Z neutrons, with each nucleon carrying the 
fraction E/A of the total kinetic energy. The relevant in¬ 
put for the initial condition is, therefore, the all-nucleon 
spectrum separated in the proton and neutron components. 
Alternatively one can use a single particle per energy bin 
and calculate particle yields at the surface for comparisons 
with full Monte Carlo methods such as E0- Fig. g 
shows the spectra of three models we typically use for cal¬ 
culations. We emphasize, that it is crucial to include the 
knee and ankle in the flux calculations with respect to the 
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Figure 2. Models of the cosmic ray nucleon spectrum and the 
neutron fraction. Gaisser-Stanev-Tilav (GST) 03 and Hillas- 
Gaisser (H3a) (H are recent 3 generation/5 mass component 
models. The proton-only broken power law model by Thunman 
et al. (TIG) G3J has been often used for calculation of the prompt 
flux in the past. The poly-gonato model ca focuses on the flux 
below and at the knee and it is not applicable at very high ener¬ 
gies. 



Figure 3. Primary model dependence of the atmospheric con¬ 
ventional + prompt neutrino flux. The model abbreviations are 
described in the caption of Fig. [2] 


range of energies accessible by current neutrino observa¬ 
tories, such as IceCube and Antares. In Fig. [3] the influ¬ 
ence of the primary model on the neutrino flux is shown. 
While at energies below the knee, where direct measure¬ 
ments of the cosmic ray flux are available, the difference 
between the models is small, there are large uncertainties 



slant depth X in [g/cm 2 ] 

Figure 4. Atmospheric density dependence on X, calculated 
using parameterizations for the US Standard Atmosphere m 
and the South Pole as implemented in CORSIKA (251 . and the 
NRLMSISE-00 model. Solid lines represent a trajectory for 
0 = 0° and dashed for 9 = 10°. 


at tens of PeV. Improving the knowledge of the spectrum 
and composition as measured by air shower experiments 
would help to disentangle these ambiguities. 


3.2 Geometry and Atmosphere 

Treating the atmosphere in planar approximation, the rela¬ 
tion between the height, slant depth, and local density can 
be taken directly from measurements. An often used ap¬ 
proach, based on the idea by Linsley IfTTI . is a parametriza- 
tion of the relation between height and mass overburden 
X v (li) (slant depth for vertical trajectory) using 5 piece- 
wise defined exponential functions, representing layers of 
the atmosphere. A higher flexibility is achieved if tabu¬ 
lated atmospheric data (e.g. from satellites) or detailed nu¬ 
merical models, such as NRLMSISE-00 (181 . are used. At 
large zenith angles also the curvature of the surface of the 
Earth has to be accounted for. Therefore we compute and 
tabulate the relation p{h atm {X)) for each provided zenith 
angle 6 and parametrization of the atmosphere. As shown 
in Fig. [4] this results in a linear smooth curve which is well 
suited for interpolation with splines. 

In Fig. [5] the ratio of the flux <1>, calculated with dif¬ 
ferent models of the atmosphere to the flux calculated us¬ 
ing the US Standard Atmosphere m is shown. Seasonal 
variations are of comparable magnitude in the CORSIKA 
parameterizations as in NRLSMSISE-00. Both models 
predict a seasonal variation of muon and neutrino rates at 
the order of +10% in agreement with what IceCube has 
observed lITil f22l. Another important feature is the non¬ 
trivial zenith angle dependence on atmospheric variations. 
A more detailed modeling of the atmosphere can be used 
to improve the experimental investigation of the prompt 
flux (231 . Furthermore, the flux at energies > PeV, which 
is dominated by prompt leptons, exhibits also a ~ 15% 
dependence on the atmosphere. 
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Muons 


Muon neutrinos 




Figure 5. Ratio of the flux calculated with different atmospheric models to the flux with US Standard atmosphere (USStd). The 
parameters and names are described in the caption of Fig. [4] The primary model is H3a and the interaction model SIBYLL-2.3 RC1. 
A vertical trajectory (9 = 0°) is represented by solid and a horizontal (6 - 90°) by dashed lines. 


4 Applications 

4.1 Calculation of the prompt flux 

In a related contribution ED, we discuss a model of 
charmed hadron production as it is implemented in the 
Monte Carlo model SIBYLL-2.3 RC1. For the calculation 
of the prompt flux using the method of this work, it is suf¬ 
ficient to take the interaction cross sections and Feynman- 
x F or the x Ulb = Secondary/E’projectiie distributions from 
Monte Carlo. Two alternative models, where it is possible 
to extract x bab distributions, are the Martin-Ryskin-Stasto 
(MRS) ll25l model and the charm model of DPMJET-II.55 
l26l . MRS considers perturbative production of charm 
quarks, based on a saturation model, and DPMJET in¬ 
cludes contributions from non-perturbative, perturbative 
and fragmentation mechanisms, similar to SIBYLL-2.3. 

The comparison in Fig. [6] shows large differences in 
shape and cross-section between the models. MRS pre¬ 
dicts a softer spectrum with a moderate forward cross- 
section. Due to a hard spectmm and the largest cross- 
section, we can expect that calculations using DPMJET 
will result in the highest fluxes. SIBYLL should produce 
more inclusive leptons compared to MRS, since the harder 
spectrum will yield charmed mesons at higher x (see dis¬ 
cussion of spectrum weighted moments in lf24l ). 

To examine the validity of charm production in DPM¬ 
JET, we compare the total cc cross-sections with ALICE 
data in Fig. [6] and the differential D ± -meson cross-section 
with forward data from LHCb in Fig. [7] The comparisons 
show that the charm model in DPMJET overestimates the 
cross-sections in the forward phase-space. DPMJET pre¬ 
dictions are disfavored by LHC data. The recent measure¬ 
ments reduce the uncertainty of charm production for pp 
collisions at PeV (Lab) energies. However, a larger frac¬ 
tion of uncertainty comes also from nuclear effects. To 
assess this uncertainty in a quantitative way, we make the 




Figure 6. (top) Feynman-.x> distributions as predicted by 
SIBYLL-2.3 RC1, the MRS model and DPMJET-II. (bottom) 
Inclusive cc cross-section in pp collsions. The ALICE measure¬ 
ment is corrected for the invisible part of the cross-section and 
extrapolated to full phase-space im 
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Figure 7. Comparison of the differential D ± cross-section with 
data, measured in pp collsions by LHCb 1281 . with SIBYLL and 
DPMJET calculations. 


assumption for the cc nuclear modification factor 


dN ce . /d p T 

p-air 1 1 1 

{N co „)&N%l&p T 


(17) 


In other words, there are no screening effects and the pro¬ 
duction of charmed quarks is a point-like process. Quan¬ 
titatively this means for the inclusive cc cross-section in 
proton-air collisions 


@~cc,p—air ~ A a j r CTc C ,p-p ~ 14.5 0~cc,p—p . (18) 

In calculations labeled SIBYLL-2.3 PL (PL for point¬ 
like) we use the charm cross sections and distributions 
of SIBYLL-2.3 RC1 as predicted for pp interactions and 
scale the yields according to Eq. [IT] 

In Fig. [ 8 ] the results of the different charm pro¬ 
duction models are compared with the Enberg-Reno- 
Sarcevic (ERS) l29l and the Thunman-Ingelman-Gondolo 
(TIG) fT5Tl calculations of the prompt muon neutrino 
flux. SIBYLL-2.3 RC1 performs similarly to ERS and it 
produces notably higher fluxes than the MRS saturation 
model. All models predict different spectral shapes. Us¬ 
ing DPMJET-II results in an order of magnitude higher 



Figure 8. Prompt muon neutrino flux calculated using models 
described in the text. To allow for direct comparison, we use the 
same primary flux model as in ERS and TIG. 


fluxes. We consider the SIBYLL-2.3 PL predictions as 
upper boundary for nuclear uncertainties of charm produc¬ 
tion. 

4.2 Partial contributions of intermediate particles 

Since fluxes of all possible intermediate mesons and 
baryons are stored in the state vector <L, we can easily 
trace back the mother particles of leptons at the surface. 
Fig. [9] is a break down of the different contributions of in¬ 
termediate particles to the conventional and prompt flux. 
We define a lepton as prompt, if its mother particle of the 
last decay has a ct < ct(K^) = 2.68 cm. To improve the 
clarity of the graph, the simple broken power-law primary 
spectrum of TIG is employed as initial state. The dom¬ 
inant contributions to conventional muons are decays of 
charged pions and kaons, while prompt muons are origi¬ 
nating from decays of charged and neutral D mesons and 
unflavored mesons, such as 77 , a> and <p. The latter reso¬ 
nances break the correlation between muon and neutrino 
fluxes at very high energies l30l as shown in the flavor 
ratios of Fig. [TO] Prompt muon neutrino and electron neu¬ 
trino fluxes are roughly equal and originate from decays 
of Z)-mesons and baryons. The fractional contribu¬ 
tion of D and AJ becomes equal at several hundreds of 
PeV. Decays of K L are responsible for the largest fraction 
of conventional muon neutrinos at energies above a few 
TeV. For electron neutrinos, other channels are important, 
such as decays of K £. At several hundreds of TeV there is 
an additional contribution from as recently discussed 
in ED. In any case, the flux from charmed particles is 
expected to be significantly higher than that due to these 
other channels. 

5 Summary and Outlook 

An efficient numerical treatment of cascade equations has 
been developed for the calculation of atmospheric lep- 
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Figure 9. Partial contribution of intermediate particles to the flux of atmospheric muons p + + p (top left), muon neutrinos (top 

right), electron neutrinos v E + v e (bottom left) and tau neutrinos v T + i> T (bottom right). The primary spectrum is TIG and the interaction 
model is SIBYLL-2.3 RC1. 



Figure 10. Flavor ratios of leptons at the surface, normalized 
to the muon neutrino flux. The calculation was performed using 
H3a primary flux and SIBYLL-2.3 RC1 for 0 = 0° (solid) and 
6 = 90° (dashed). 


ton fluxes at very high energy, with particular attention 
put on the transition from conventional to prompt produc¬ 
tion processes. As a first application we have shown cal¬ 
culations to illustrate the importance of the primary flux 
parametrization, the model of the atmosphere, the role of 


short-lived particles, and the model of charmed hadron 
production. More details about the numerical solution and 
the code will be published elsewhere. 
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